Summarise stomach content data

Author

Max Lindmark

Published

December 11, 2023

Load packages

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "RCurl", "janitor", "terra", "viridis", "devtools", "ncdf4", "chron", "raster")

if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)

}

invisible(lapply(pkgs, library, character.only = T))

# Source code for map plots
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")

home <- here::here()

Read data

d <- read_csv(paste0(home, "/data/clean/full_stomach_data.csv"))

Summarize and organize data

We want 1 row = 1 predator and the total weight for each present prey type

# Calculate total weight of prey by predator ID and prey species (i.e., across prey sizes). First create wide data frame so that I can sum easily across prey groups (columns)
d_wide <- d %>% 
  group_by(pred_id, prey_latin_name) %>% 
  summarise(tot_prey_weight_g = sum(prey_weight_g)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = prey_latin_name, values_from = tot_prey_weight_g) %>% 
  mutate_all(~ifelse(is.na(.), 0, .)) %>% 
  clean_names()
  
# Now make some calculations and aggregate to some taxonomic level. Since all columns are assigned to some higher level group (or the same group), the sum of these is the total stomach content. Note that I have one group for unidentified clupeids, but also sprat and herring. So if I want the total of some aggregated group, then I need to add all the sub-groups.

sort(colnames(d_wide %>% dplyr::select(-pred_id)))
  [1] "aglae"                    "agonus_cataphractus"     
  [3] "algae"                    "ammodytidae"             
  [5] "amphipoda"                "aphia_minuta"            
  [7] "bivalvia"                 "bylgides_sarsi"          
  [9] "carbon"                   "carbon_2"                
 [11] "caridea"                  "clupea_harengus"         
 [13] "clupea_harengus_2"        "clupeidae"               
 [15] "clupeidae_2"              "copepoda"                
 [17] "crangon"                  "crangon_crangon"         
 [19] "crustacea"                "cumacea"                 
 [21] "decapoda"                 "diastylis_rathkei"       
 [23] "digestive_tract"          "enchelyopus_cimbrius"    
 [25] "engraulis_encrasicolus"   "gadus_morhua"            
 [27] "gammarus_sp"              "gasterosteidae"          
 [29] "gasterosteus_aculeatus"   "gastrosacus"             
 [31] "gobiidae"                 "gobiidae_2"              
 [33] "gobius_niger"             "halicryptus"             
 [35] "halicryptus_spinolusus"   "halicryptus_spinulosus"  
 [37] "halicryptus_spinulosus_2" "hediste_diversicolor"    
 [39] "hydrobia_sp"              "hyperia_galba"           
 [41] "idotea_balthica"          "idotea_sp"               
 [43] "limecola_balthica"        "limecola_balthica_2"     
 [45] "litter_plastics"          "macoma_balthica"         
 [47] "monoporeia_affinis"       "mucus"                   
 [49] "mucus_2"                  "mya_arenaria"            
 [51] "myoxocephalus_scorpius"   "mysida"                  
 [53] "mysidae"                  "mysis_mixta"             
 [55] "mytilus_edulis"           "mytilus_sp"              
 [57] "mytilus_sp_2"             "na"                      
 [59] "neogobius_melanostomus"   "neomysis_integer"        
 [61] "nephtys_ciliata"          "palaemon_elegans"        
 [63] "palaemon_sp"              "pectinaria_sp"           
 [65] "phyllodocida"             "pisces"                  
 [67] "pisces_2"                 "plastic"                 
 [69] "plastics"                 "platichthys_flesus"      
 [71] "pleuronectidae"           "polychaeta"              
 [73] "pontoporeia_femorata"     "pontoporeia_femotara"    
 [75] "pontoporeidae"            "pontoporeiidae"          
 [77] "prapulida"                "praunus_flexuosus"       
 [79] "priapulida"               "priapulida_2"            
 [81] "priapulus"                "priapulus_caudatus"      
 [83] "priapulus_caudatus_2"     "pungitius_pungitius"     
 [85] "remains"                  "remains_2"               
 [87] "saduria_entomon"          "sand"                    
 [89] "sand_2"                   "scales"                  
 [91] "scales_2"                 "scoloplos_armiger"       
 [93] "scomber_scombrus"         "sprattus_sprattus"       
 [95] "sprattus_sprattus_2"      "stone"                   
 [97] "stone_2"                  "waste"                   
 [99] "wood"                     "wood_2"                  
[101] "zoarces_viviparus"       
d_wide2 <- d_wide %>% 
  mutate(amphipoda_tot = gammarus_sp + monoporeia_affinis + amphipoda + hyperia_galba,
         bivalvia_tot =  bivalvia + mytilus_sp + mytilus_sp_2 + mya_arenaria + macoma_balthica + 
             mytilus_edulis + limecola_balthica + limecola_balthica_2,
         clupeidae_tot = clupeidae + clupeidae_2,
         clupea_harengus_tot = clupea_harengus + clupea_harengus_2,
         gadus_morhua_tot = gadus_morhua,
         gobiidae_tot = gobiidae + gobiidae_2 + gobius_niger + neogobius_melanostomus + aphia_minuta,
         mysidae_tot = mysidae + neomysis_integer + mysis_mixta + mysida + gastrosacus,
         non_bio_tot = stone + stone_2 + plastic + plastics + sand + wood + carbon + carbon_2 + wood_2 + 
             litter_plastics + sand_2,
         other_crustacea_tot = pontoporeia_femorata + pontoporeia_femotara + crangon + 
             crangon_crangon + idotea_balthica + cumacea + idotea_sp +
             praunus_flexuosus + crustacea + diastylis_rathkei + palaemon_sp + palaemon_elegans + caridea +
             copepoda + pontoporeiidae + decapoda + pontoporeidae, 
         other_tot = halicryptus_spinulosus + halicryptus_spinulosus_2 + priapulus + priapulus_caudatus + priapulus_caudatus_2 + algae + aglae + 
             waste + remains + remains_2 + hydrobia_sp + 
             priapulida + halicryptus + digestive_tract + mucus + mucus_2 + 
             halicryptus_spinolusus + priapulida_2 + prapulida,
         other_pisces_tot = pisces + pisces_2 + enchelyopus_cimbrius + engraulis_encrasicolus +
             gasterosteus_aculeatus + scales + scales_2 + scomber_scombrus +
             pungitius_pungitius + zoarces_viviparus + ammodytidae +
             pleuronectidae + gasterosteidae + agonus_cataphractus + myoxocephalus_scorpius,
         platichthys_flesus_tot = platichthys_flesus,
         polychaeta_tot = bylgides_sarsi + scoloplos_armiger + hediste_diversicolor +
             phyllodocida + polychaeta + pectinaria_sp + nephtys_ciliata,
         saduria_entomon_tot = saduria_entomon,
         sprattus_sprattus_tot = sprattus_sprattus + sprattus_sprattus_2
         )

# Check it's correct! The sum of all groups should be the sum of all categories
colnames(d_wide2 %>% dplyr::select(-pred_id))
  [1] "diastylis_rathkei"        "halicryptus_spinulosus"  
  [3] "bylgides_sarsi"           "na"                      
  [5] "clupeidae"                "sprattus_sprattus"       
  [7] "mysis_mixta"              "stone"                   
  [9] "crangon_crangon"          "gammarus_sp"             
 [11] "priapulida"               "priapulus_caudatus"      
 [13] "gasterosteus_aculeatus"   "pisces"                  
 [15] "saduria_entomon"          "neomysis_integer"        
 [17] "clupea_harengus"          "monoporeia_affinis"      
 [19] "waste"                    "scales"                  
 [21] "gobiidae"                 "pontoporeia_femorata"    
 [23] "remains"                  "crustacea"               
 [25] "mysidae"                  "limecola_balthica"       
 [27] "bivalvia"                 "halicryptus"             
 [29] "zoarces_viviparus"        "platichthys_flesus"      
 [31] "digestive_tract"          "wood"                    
 [33] "phyllodocida"             "polychaeta"              
 [35] "sand"                     "mytilus_sp"              
 [37] "amphipoda"                "limecola_balthica_2"     
 [39] "algae"                    "pungitius_pungitius"     
 [41] "scoloplos_armiger"        "gadus_morhua"            
 [43] "priapulida_2"             "idotea_sp"               
 [45] "enchelyopus_cimbrius"     "pleuronectidae"          
 [47] "cumacea"                  "plastic"                 
 [49] "stone_2"                  "crangon"                 
 [51] "sprattus_sprattus_2"      "aglae"                   
 [53] "macoma_balthica"          "carbon"                  
 [55] "gasterosteidae"           "mysida"                  
 [57] "gobiidae_2"               "gobius_niger"            
 [59] "palaemon_sp"              "mytilus_sp_2"            
 [61] "scales_2"                 "ammodytidae"             
 [63] "pectinaria_sp"            "sand_2"                  
 [65] "pontoporeiidae"           "mucus"                   
 [67] "pontoporeia_femotara"     "remains_2"               
 [69] "mucus_2"                  "priapulus"               
 [71] "carbon_2"                 "wood_2"                  
 [73] "halicryptus_spinulosus_2" "pisces_2"                
 [75] "mya_arenaria"             "gastrosacus"             
 [77] "nephtys_ciliata"          "litter_plastics"         
 [79] "clupeidae_2"              "pontoporeidae"           
 [81] "decapoda"                 "praunus_flexuosus"       
 [83] "neogobius_melanostomus"   "plastics"                
 [85] "agonus_cataphractus"      "clupea_harengus_2"       
 [87] "copepoda"                 "halicryptus_spinolusus"  
 [89] "prapulida"                "mytilus_edulis"          
 [91] "hydrobia_sp"              "myoxocephalus_scorpius"  
 [93] "idotea_balthica"          "caridea"                 
 [95] "aphia_minuta"             "hediste_diversicolor"    
 [97] "palaemon_elegans"         "priapulus_caudatus_2"    
 [99] "scomber_scombrus"         "hyperia_galba"           
[101] "engraulis_encrasicolus"   "amphipoda_tot"           
[103] "bivalvia_tot"             "clupeidae_tot"           
[105] "clupea_harengus_tot"      "gadus_morhua_tot"        
[107] "gobiidae_tot"             "mysidae_tot"             
[109] "non_bio_tot"              "other_crustacea_tot"     
[111] "other_tot"                "other_pisces_tot"        
[113] "platichthys_flesus_tot"   "polychaeta_tot"          
[115] "saduria_entomon_tot"      "sprattus_sprattus_tot"   
t <- d_wide2 %>%
  dplyr::select(-pred_id) %>% 
  mutate(all_cat_sum = rowSums(.[1:100]),
         group_sum = rowSums(.[101:115])) %>% 
  dplyr::select(all_cat_sum, group_sum) %>%
  filter(!all_cat_sum == group_sum)  

# Must be some rounding error...
d_wide2 %>%
  dplyr::select(-pred_id) %>% 
  mutate(all_cat_sum = rowSums(.[1:100]),
         group_sum = rowSums(.[101:115])) %>% 
  dplyr::select(all_cat_sum, group_sum) %>%
  ggplot(aes(all_cat_sum, group_sum)) +
  geom_point() +
  geom_abline(color = "red")

# Select only columns aggregated columns (ending with _tot) (all columns (prey) are represented there)
colnames(d_wide2)
  [1] "pred_id"                  "diastylis_rathkei"       
  [3] "halicryptus_spinulosus"   "bylgides_sarsi"          
  [5] "na"                       "clupeidae"               
  [7] "sprattus_sprattus"        "mysis_mixta"             
  [9] "stone"                    "crangon_crangon"         
 [11] "gammarus_sp"              "priapulida"              
 [13] "priapulus_caudatus"       "gasterosteus_aculeatus"  
 [15] "pisces"                   "saduria_entomon"         
 [17] "neomysis_integer"         "clupea_harengus"         
 [19] "monoporeia_affinis"       "waste"                   
 [21] "scales"                   "gobiidae"                
 [23] "pontoporeia_femorata"     "remains"                 
 [25] "crustacea"                "mysidae"                 
 [27] "limecola_balthica"        "bivalvia"                
 [29] "halicryptus"              "zoarces_viviparus"       
 [31] "platichthys_flesus"       "digestive_tract"         
 [33] "wood"                     "phyllodocida"            
 [35] "polychaeta"               "sand"                    
 [37] "mytilus_sp"               "amphipoda"               
 [39] "limecola_balthica_2"      "algae"                   
 [41] "pungitius_pungitius"      "scoloplos_armiger"       
 [43] "gadus_morhua"             "priapulida_2"            
 [45] "idotea_sp"                "enchelyopus_cimbrius"    
 [47] "pleuronectidae"           "cumacea"                 
 [49] "plastic"                  "stone_2"                 
 [51] "crangon"                  "sprattus_sprattus_2"     
 [53] "aglae"                    "macoma_balthica"         
 [55] "carbon"                   "gasterosteidae"          
 [57] "mysida"                   "gobiidae_2"              
 [59] "gobius_niger"             "palaemon_sp"             
 [61] "mytilus_sp_2"             "scales_2"                
 [63] "ammodytidae"              "pectinaria_sp"           
 [65] "sand_2"                   "pontoporeiidae"          
 [67] "mucus"                    "pontoporeia_femotara"    
 [69] "remains_2"                "mucus_2"                 
 [71] "priapulus"                "carbon_2"                
 [73] "wood_2"                   "halicryptus_spinulosus_2"
 [75] "pisces_2"                 "mya_arenaria"            
 [77] "gastrosacus"              "nephtys_ciliata"         
 [79] "litter_plastics"          "clupeidae_2"             
 [81] "pontoporeidae"            "decapoda"                
 [83] "praunus_flexuosus"        "neogobius_melanostomus"  
 [85] "plastics"                 "agonus_cataphractus"     
 [87] "clupea_harengus_2"        "copepoda"                
 [89] "halicryptus_spinolusus"   "prapulida"               
 [91] "mytilus_edulis"           "hydrobia_sp"             
 [93] "myoxocephalus_scorpius"   "idotea_balthica"         
 [95] "caridea"                  "aphia_minuta"            
 [97] "hediste_diversicolor"     "palaemon_elegans"        
 [99] "priapulus_caudatus_2"     "scomber_scombrus"        
[101] "hyperia_galba"            "engraulis_encrasicolus"  
[103] "amphipoda_tot"            "bivalvia_tot"            
[105] "clupeidae_tot"            "clupea_harengus_tot"     
[107] "gadus_morhua_tot"         "gobiidae_tot"            
[109] "mysidae_tot"              "non_bio_tot"             
[111] "other_crustacea_tot"      "other_tot"               
[113] "other_pisces_tot"         "platichthys_flesus_tot"  
[115] "polychaeta_tot"           "saduria_entomon_tot"     
[117] "sprattus_sprattus_tot"   
d_wide3 <- d_wide2 %>%
  dplyr::select(pred_id, ends_with("_tot"))

# Add back in other information about the predator ID
d_sel <- d %>%
  dplyr::select(predator_latin_name, species, pred_weight_g, pred_length_cm,
                year, quarter, month, ices_rect, subdiv, haul_id, smhi_serial_no,
                X, Y, lat, lon, pred_id, fle_kg_km2, mcod_kg_km2, scod_kg_km2, date,
                bottom_depth
                ) %>% 
  distinct(pred_id, .keep_all = TRUE)

d_wide3 <- left_join(d_wide3, d_sel)

# mutate(size_group = ifelse(lengthcl >= 100 & lengthcl < 250, "small", NA),
#        size_group = ifelse(lengthcl >= 250 & lengthcl < 500, "medium", size_group)) %>% 

d_wide3 <- d_wide3 %>%
  mutate(group = NA,
         group = ifelse(species == "Flounder", "flounder", group),
         group = ifelse(species == "Cod" & pred_length_cm >= 10 & pred_length_cm <= 25, "small cod", group),
         group = ifelse(species == "Cod" & pred_length_cm > 25 & pred_length_cm <= 50, "large cod", group))

Add covariates

Depth

dep_raster <- terra::rast(paste0(home, "/data/Mean depth natural colour (with land).nc"))
class(dep_raster)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
crs(dep_raster, proj = TRUE)
[1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)

d_wide3$depth <- terra::extract(dep_raster, d_wide3 %>% dplyr::select(lon, lat))$elevation

ggplot(d_wide3, aes(lon, lat, color = depth*-1)) + 
  geom_point()

d_wide3$depth <- d_wide3$depth*-1

d_wide3 <- d_wide3 |> drop_na(depth)
drop_na: no rows removed
plot_map + 
  geom_point(data = d_wide3, aes(X*1000, Y*1000, color = depth)) + 
  scale_color_viridis()

ggplot(d_wide3, aes(depth, bottom_depth)) + 
  geom_point() + 
  geom_abline()

# So strange. What does it look like in the normal data?
t <- d_wide3 %>% 
  filter(bottom_depth > 250)
filter: removed 5,396 rows (56%), 4,253 rows remaining
ggplot(t, aes(depth, bottom_depth/10)) + 
  geom_point() + 
  geom_abline()

strange_dept <- d %>% filter(haul_id %in% c(t$haul_id))
filter: removed 49,130 rows (54%), 42,145 rows remaining
strange_dept %>% 
  summary()
   pred_id            data_id          gall_bladder       stomach_state     
 Length:42145       Length:42145       Length:42145       Length:42145      
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 predator_code      prey_weight_g       stage_digestion   comment         
 Length:42145       Min.   :  0.00000   Min.   :0.000   Length:42145      
 Class :character   1st Qu.:  0.01069   1st Qu.:1.000   Class :character  
 Mode  :character   Median :  0.02562   Median :1.000   Mode  :character  
                    Mean   :  0.31695   Mean   :1.063                     
                    3rd Qu.:  0.06588   3rd Qu.:1.000                     
                    Max.   :144.34000   Max.   :2.000                     
                                        NA's   :1118                      
 prey_length_cm   prey_latin_name    length_code           n_empty     
 Min.   :  4.00   Length:42145       Length:42145       Min.   :0.000  
 1st Qu.: 20.00   Class :character   Class :character   1st Qu.:0.000  
 Median : 30.00   Mode  :character   Mode  :character   Median :0.000  
 Mean   : 45.83                                         Mean   :0.031  
 3rd Qu.: 48.00                                         3rd Qu.:0.000  
 Max.   :300.00                                         Max.   :1.000  
 NA's   :39172                                          NA's   :10698  
  n_with_food     n_regurgitated    haul_id            fle_kg_km2     
 Min.   :0.0000   Min.   :0.000   Length:42145       Min.   :    0.0  
 1st Qu.:1.0000   1st Qu.:0.000   Class :character   1st Qu.:  540.5  
 Median :1.0000   Median :0.000   Mode  :character   Median : 1137.2  
 Mean   :0.9761   Mean   :0.003                      Mean   : 1740.3  
 3rd Qu.:1.0000   3rd Qu.:0.000                      3rd Qu.: 1961.9  
 Max.   :1.0000   Max.   :1.000                      Max.   :18668.7  
 NA's   :313      NA's   :10911                                       
  mcod_kg_km2        scod_kg_km2           year         quarter     
 Min.   :     0.0   Min.   :    0.0   Min.   :2015   Min.   :1.000  
 1st Qu.:   149.6   1st Qu.:  102.0   1st Qu.:2017   1st Qu.:1.000  
 Median :  1287.8   Median :  462.2   Median :2021   Median :4.000  
 Mean   :  4724.4   Mean   : 1771.7   Mean   :2020   Mean   :3.184  
 3rd Qu.:  4007.8   3rd Qu.: 2012.0   3rd Qu.:2022   3rd Qu.:4.000  
 Max.   :171072.0   Max.   :41057.3   Max.   :2022   Max.   :4.000  
                                                                    
      lat             lon         bottom_depth      species         
 Min.   :55.17   Min.   :13.30   Min.   : 333.0   Length:42145      
 1st Qu.:55.67   1st Qu.:14.47   1st Qu.: 418.0   Class :character  
 Median :55.82   Median :15.80   Median : 473.0   Mode  :character  
 Mean   :56.24   Mean   :16.31   Mean   : 490.7                     
 3rd Qu.:57.15   3rd Qu.:18.40   3rd Qu.: 552.0                     
 Max.   :58.45   Max.   :19.52   Max.   :1216.0                     
                                                                    
 pred_length_cm  pred_weight_g         size       subdiv     
 Min.   : 4.00   Min.   :   1.0   Min.   :9   Min.   :24.00  
 1st Qu.:23.00   1st Qu.: 120.0   1st Qu.:9   1st Qu.:25.00  
 Median :26.00   Median : 185.0   Median :9   Median :25.00  
 Mean   :26.66   Mean   : 224.9   Mean   :9   Mean   :25.74  
 3rd Qu.:31.00   3rd Qu.: 292.0   3rd Qu.:9   3rd Qu.:27.00  
 Max.   :77.00   Max.   :4798.0   Max.   :9   Max.   :28.00  
                                                             
      date                month              X               Y       
 Min.   :2015-11-20   Min.   : 2.000   Min.   :391.7   Min.   :6115  
 1st Qu.:2017-02-26   1st Qu.: 3.000   1st Qu.:466.5   1st Qu.:6169  
 Median :2021-11-21   Median :11.000   Median :550.2   Median :6186  
 Mean   :2020-05-24   Mean   : 8.588   Mean   :578.5   Mean   :6236  
 3rd Qu.:2022-02-26   3rd Qu.:11.000   3rd Qu.:710.3   3rd Qu.:6341  
 Max.   :2022-11-27   Max.   :11.000   Max.   :768.0   Max.   :6487  
                                                                     
 smhi_serial_no   prey_number_type   prey_weight_type    ices_rect        
 Min.   :   1.0   Length:42145       Length:42145       Length:42145      
 1st Qu.:  17.0   Class :character   Class :character   Class :character  
 Median : 201.0   Mode  :character   Mode  :character   Mode  :character  
 Mean   : 442.2                                                           
 3rd Qu.: 969.0                                                           
 Max.   :1014.0                                                           
 NA's   :13421                                                            
 predator_latin_name   rel_weight       
 Length:42145        Min.   :0.000e+00  
 Class :character    1st Qu.:6.357e-05  
 Mode  :character    Median :1.657e-04  
                     Mean   :9.743e-04  
                     3rd Qu.:4.557e-04  
                     Max.   :8.300e-02  
                                        
unique(strange_dept$quarter)
[1] 4 1
unique(strange_dept$year)
[1] 2015 2016 2017 2019 2020 2021 2022
unique(strange_dept$month)
[1] 11  2  3
unique(strange_dept$subdiv)
[1] 25 26 28 24 27

Add in environmental variables

dat <- d_wide3 %>% 
  dplyr::select(haul_id, X, Y, lon, lat, year, month, quarter) %>% 
  distinct(haul_id, .keep_all = TRUE) %>% 
  mutate(month = ifelse(quarter == 1, 2, 11), # most common months
         month_year = paste(month, year, sep = "_"))
distinct: removed 9,295 rows (96%), 354 rows remaining
mutate: changed 78 values (22%) of 'month' (0 new NA)
        new variable 'month_year' (character) with 14 unique values and 0% NA
dat %>% head()
# A tibble: 6 × 9
  haul_id       X     Y   lon   lat  year month quarter month_year
  <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <chr>     
1 2015_4_2   475. 6165.  14.6  55.6  2015    11       4 11_2015   
2 2015_4_4   469. 6171.  14.5  55.7  2015    11       4 11_2015   
3 2015_4_6   460. 6173.  14.4  55.7  2015    11       4 11_2015   
4 2015_4_8   647. 6216.  17.4  56.1  2015    11       4 11_2015   
5 2015_4_9   710. 6247.  18.4  56.3  2015    11       4 11_2015   
6 2015_4_12  722. 6255.  18.6  56.4  2015    11       4 11_2015   
# For temperature and oxygen, we need to use two different rasters because we don't have the updated data in the re-analyisis
# TODO: check this later
dat_1 <- dat %>% filter(year < 2022)
filter: removed 49 rows (14%), 305 rows remaining
dat_2 <- dat %>% filter(year >= 2022)
filter: removed 305 rows (86%), 49 rows remaining

Saduria

saduria <- terra::rast(paste0(home, "/data/saduria-tif/FWBiomassm_raster_19812019presHighweightcor_no0_newZi.tif"))

WGS84 <- "+proj=longlat +datum=WGS84"

saduria_latlon <- terra::project(saduria, WGS84)

density_saduria <- terra::extract(saduria_latlon, dat %>% dplyr::select(lon, lat))

dat$density_saduria <- density_saduria$FWBiomassm_raster_19812019presHighweightcor_no0_newZi

plot_map + 
  geom_point(data = dat, aes(X*1000, Y*1000, color = density_saduria)) +
  scale_color_viridis(trans = "sqrt", name = "Saduria biomass density")

Oxygen

# Downloaded from here: https://data.marine.copernicus.eu/products
# Extract raster points: https://gisday.wordpress.com/2014/03/24/extract-raster-values-from-points-using-r/comment-page-1/
# https://rpubs.com/boyerag/297592
# https://pjbartlein.github.io/REarthSysSci/netCDF.html#get-a-variable
# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_bgc_my_P1D-m_1702039771395.nc"))

print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/NEMO/cmems_mod_bal_bgc_my_P1D-m_1702039771395.nc (NC_FORMAT_CLASSIC):

     1 variables (excluding dimension variables):
        float o2b[lon,lat,time]   
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            long_name: Dissolved_Oxygen_Concentration
            units: mmol m-3
            _FillValue: -999
            missing_value: -999
            unit_long: millimole O2 per unit volume
            interval_write: 1 d
            interval_operation: 90 s
            online_operation: average
            _ChunkSizes: 1
             _ChunkSizes: 774
             _ChunkSizes: 763

     3 dimensions:
        time  Size:2557 
            standard_name: time
            long_name: time
            units: days since 1900-01-01 00:00:00
            calendar: standard
            axis: T
            _ChunkSizes: 512
            _CoordinateAxisType: Time
            valid_min: 42003.5
            valid_max: 44559.5
        lat  Size:289 
            standard_name: latitude
            long_name: Latitude of each location
            units: degrees_north
            axis: Y
            _CoordinateAxisType: Lat
            valid_min: 53.2582931518555
            valid_max: 58.0582122802734
        lon  Size:336 
            standard_name: longitude
            long_name: Longitude of each location
            units: degrees_east
            axis: X
            _CoordinateAxisType: Lon
            valid_min: 12.5414085388184
            valid_max: 21.8466091156006

    22 global attributes:
        CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
        Conventions: CF-1.0
        history: Tue Mar 14 17:24:24 2023: cdo -z zip_1 --sortname copy /net/isilon/ifs/arch/home/iri/BALMFC/Nemo4.0/BALMFCMYP2022_univariateSST//BAL-ERGOM_BGC-DailyMeans-20211231.nc tmp_20211231.nc
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        institution: Baltic MFC, PU Danish Meteorological Institute
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        title: CMEMS ERGOM daily integrated model fields
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        compression: yes
        stop_date: 2021-12-31 12:00:00
        creation_date: 2023-03-14 17:24:23
        FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
        FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
        FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
        FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
        FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
        start_date: 2021-12-31 12:00:00
        CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
        _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 12.54141 12.56919 12.59696 12.62474 12.65252 12.68029
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.25829 53.27496 53.29163 53.30829 53.32496 53.34163
# Get time
time <- ncvar_get(ncin,"time")
time
   [1] 42003.5 42004.5 42005.5 42006.5 42007.5 42008.5 42009.5 42010.5 42011.5
  [10] 42012.5 42013.5 42014.5 42015.5 42016.5 42017.5 42018.5 42019.5 42020.5
  [19] 42021.5 42022.5 42023.5 42024.5 42025.5 42026.5 42027.5 42028.5 42029.5
  [28] 42030.5 42031.5 42032.5 42033.5 42034.5 42035.5 42036.5 42037.5 42038.5
  [37] 42039.5 42040.5 42041.5 42042.5 42043.5 42044.5 42045.5 42046.5 42047.5
  [46] 42048.5 42049.5 42050.5 42051.5 42052.5 42053.5 42054.5 42055.5 42056.5
  [55] 42057.5 42058.5 42059.5 42060.5 42061.5 42062.5 42063.5 42064.5 42065.5
  [64] 42066.5 42067.5 42068.5 42069.5 42070.5 42071.5 42072.5 42073.5 42074.5
  [73] 42075.5 42076.5 42077.5 42078.5 42079.5 42080.5 42081.5 42082.5 42083.5
  [82] 42084.5 42085.5 42086.5 42087.5 42088.5 42089.5 42090.5 42091.5 42092.5
  [91] 42093.5 42094.5 42095.5 42096.5 42097.5 42098.5 42099.5 42100.5 42101.5
 [100] 42102.5 42103.5 42104.5 42105.5 42106.5 42107.5 42108.5 42109.5 42110.5
 [109] 42111.5 42112.5 42113.5 42114.5 42115.5 42116.5 42117.5 42118.5 42119.5
 [118] 42120.5 42121.5 42122.5 42123.5 42124.5 42125.5 42126.5 42127.5 42128.5
 [127] 42129.5 42130.5 42131.5 42132.5 42133.5 42134.5 42135.5 42136.5 42137.5
 [136] 42138.5 42139.5 42140.5 42141.5 42142.5 42143.5 42144.5 42145.5 42146.5
 [145] 42147.5 42148.5 42149.5 42150.5 42151.5 42152.5 42153.5 42154.5 42155.5
 [154] 42156.5 42157.5 42158.5 42159.5 42160.5 42161.5 42162.5 42163.5 42164.5
 [163] 42165.5 42166.5 42167.5 42168.5 42169.5 42170.5 42171.5 42172.5 42173.5
 [172] 42174.5 42175.5 42176.5 42177.5 42178.5 42179.5 42180.5 42181.5 42182.5
 [181] 42183.5 42184.5 42185.5 42186.5 42187.5 42188.5 42189.5 42190.5 42191.5
 [190] 42192.5 42193.5 42194.5 42195.5 42196.5 42197.5 42198.5 42199.5 42200.5
 [199] 42201.5 42202.5 42203.5 42204.5 42205.5 42206.5 42207.5 42208.5 42209.5
 [208] 42210.5 42211.5 42212.5 42213.5 42214.5 42215.5 42216.5 42217.5 42218.5
 [217] 42219.5 42220.5 42221.5 42222.5 42223.5 42224.5 42225.5 42226.5 42227.5
 [226] 42228.5 42229.5 42230.5 42231.5 42232.5 42233.5 42234.5 42235.5 42236.5
 [235] 42237.5 42238.5 42239.5 42240.5 42241.5 42242.5 42243.5 42244.5 42245.5
 [244] 42246.5 42247.5 42248.5 42249.5 42250.5 42251.5 42252.5 42253.5 42254.5
 [253] 42255.5 42256.5 42257.5 42258.5 42259.5 42260.5 42261.5 42262.5 42263.5
 [262] 42264.5 42265.5 42266.5 42267.5 42268.5 42269.5 42270.5 42271.5 42272.5
 [271] 42273.5 42274.5 42275.5 42276.5 42277.5 42278.5 42279.5 42280.5 42281.5
 [280] 42282.5 42283.5 42284.5 42285.5 42286.5 42287.5 42288.5 42289.5 42290.5
 [289] 42291.5 42292.5 42293.5 42294.5 42295.5 42296.5 42297.5 42298.5 42299.5
 [298] 42300.5 42301.5 42302.5 42303.5 42304.5 42305.5 42306.5 42307.5 42308.5
 [307] 42309.5 42310.5 42311.5 42312.5 42313.5 42314.5 42315.5 42316.5 42317.5
 [316] 42318.5 42319.5 42320.5 42321.5 42322.5 42323.5 42324.5 42325.5 42326.5
 [325] 42327.5 42328.5 42329.5 42330.5 42331.5 42332.5 42333.5 42334.5 42335.5
 [334] 42336.5 42337.5 42338.5 42339.5 42340.5 42341.5 42342.5 42343.5 42344.5
 [343] 42345.5 42346.5 42347.5 42348.5 42349.5 42350.5 42351.5 42352.5 42353.5
 [352] 42354.5 42355.5 42356.5 42357.5 42358.5 42359.5 42360.5 42361.5 42362.5
 [361] 42363.5 42364.5 42365.5 42366.5 42367.5 42368.5 42369.5 42370.5 42371.5
 [370] 42372.5 42373.5 42374.5 42375.5 42376.5 42377.5 42378.5 42379.5 42380.5
 [379] 42381.5 42382.5 42383.5 42384.5 42385.5 42386.5 42387.5 42388.5 42389.5
 [388] 42390.5 42391.5 42392.5 42393.5 42394.5 42395.5 42396.5 42397.5 42398.5
 [397] 42399.5 42400.5 42401.5 42402.5 42403.5 42404.5 42405.5 42406.5 42407.5
 [406] 42408.5 42409.5 42410.5 42411.5 42412.5 42413.5 42414.5 42415.5 42416.5
 [415] 42417.5 42418.5 42419.5 42420.5 42421.5 42422.5 42423.5 42424.5 42425.5
 [424] 42426.5 42427.5 42428.5 42429.5 42430.5 42431.5 42432.5 42433.5 42434.5
 [433] 42435.5 42436.5 42437.5 42438.5 42439.5 42440.5 42441.5 42442.5 42443.5
 [442] 42444.5 42445.5 42446.5 42447.5 42448.5 42449.5 42450.5 42451.5 42452.5
 [451] 42453.5 42454.5 42455.5 42456.5 42457.5 42458.5 42459.5 42460.5 42461.5
 [460] 42462.5 42463.5 42464.5 42465.5 42466.5 42467.5 42468.5 42469.5 42470.5
 [469] 42471.5 42472.5 42473.5 42474.5 42475.5 42476.5 42477.5 42478.5 42479.5
 [478] 42480.5 42481.5 42482.5 42483.5 42484.5 42485.5 42486.5 42487.5 42488.5
 [487] 42489.5 42490.5 42491.5 42492.5 42493.5 42494.5 42495.5 42496.5 42497.5
 [496] 42498.5 42499.5 42500.5 42501.5 42502.5 42503.5 42504.5 42505.5 42506.5
 [505] 42507.5 42508.5 42509.5 42510.5 42511.5 42512.5 42513.5 42514.5 42515.5
 [514] 42516.5 42517.5 42518.5 42519.5 42520.5 42521.5 42522.5 42523.5 42524.5
 [523] 42525.5 42526.5 42527.5 42528.5 42529.5 42530.5 42531.5 42532.5 42533.5
 [532] 42534.5 42535.5 42536.5 42537.5 42538.5 42539.5 42540.5 42541.5 42542.5
 [541] 42543.5 42544.5 42545.5 42546.5 42547.5 42548.5 42549.5 42550.5 42551.5
 [550] 42552.5 42553.5 42554.5 42555.5 42556.5 42557.5 42558.5 42559.5 42560.5
 [559] 42561.5 42562.5 42563.5 42564.5 42565.5 42566.5 42567.5 42568.5 42569.5
 [568] 42570.5 42571.5 42572.5 42573.5 42574.5 42575.5 42576.5 42577.5 42578.5
 [577] 42579.5 42580.5 42581.5 42582.5 42583.5 42584.5 42585.5 42586.5 42587.5
 [586] 42588.5 42589.5 42590.5 42591.5 42592.5 42593.5 42594.5 42595.5 42596.5
 [595] 42597.5 42598.5 42599.5 42600.5 42601.5 42602.5 42603.5 42604.5 42605.5
 [604] 42606.5 42607.5 42608.5 42609.5 42610.5 42611.5 42612.5 42613.5 42614.5
 [613] 42615.5 42616.5 42617.5 42618.5 42619.5 42620.5 42621.5 42622.5 42623.5
 [622] 42624.5 42625.5 42626.5 42627.5 42628.5 42629.5 42630.5 42631.5 42632.5
 [631] 42633.5 42634.5 42635.5 42636.5 42637.5 42638.5 42639.5 42640.5 42641.5
 [640] 42642.5 42643.5 42644.5 42645.5 42646.5 42647.5 42648.5 42649.5 42650.5
 [649] 42651.5 42652.5 42653.5 42654.5 42655.5 42656.5 42657.5 42658.5 42659.5
 [658] 42660.5 42661.5 42662.5 42663.5 42664.5 42665.5 42666.5 42667.5 42668.5
 [667] 42669.5 42670.5 42671.5 42672.5 42673.5 42674.5 42675.5 42676.5 42677.5
 [676] 42678.5 42679.5 42680.5 42681.5 42682.5 42683.5 42684.5 42685.5 42686.5
 [685] 42687.5 42688.5 42689.5 42690.5 42691.5 42692.5 42693.5 42694.5 42695.5
 [694] 42696.5 42697.5 42698.5 42699.5 42700.5 42701.5 42702.5 42703.5 42704.5
 [703] 42705.5 42706.5 42707.5 42708.5 42709.5 42710.5 42711.5 42712.5 42713.5
 [712] 42714.5 42715.5 42716.5 42717.5 42718.5 42719.5 42720.5 42721.5 42722.5
 [721] 42723.5 42724.5 42725.5 42726.5 42727.5 42728.5 42729.5 42730.5 42731.5
 [730] 42732.5 42733.5 42734.5 42735.5 42736.5 42737.5 42738.5 42739.5 42740.5
 [739] 42741.5 42742.5 42743.5 42744.5 42745.5 42746.5 42747.5 42748.5 42749.5
 [748] 42750.5 42751.5 42752.5 42753.5 42754.5 42755.5 42756.5 42757.5 42758.5
 [757] 42759.5 42760.5 42761.5 42762.5 42763.5 42764.5 42765.5 42766.5 42767.5
 [766] 42768.5 42769.5 42770.5 42771.5 42772.5 42773.5 42774.5 42775.5 42776.5
 [775] 42777.5 42778.5 42779.5 42780.5 42781.5 42782.5 42783.5 42784.5 42785.5
 [784] 42786.5 42787.5 42788.5 42789.5 42790.5 42791.5 42792.5 42793.5 42794.5
 [793] 42795.5 42796.5 42797.5 42798.5 42799.5 42800.5 42801.5 42802.5 42803.5
 [802] 42804.5 42805.5 42806.5 42807.5 42808.5 42809.5 42810.5 42811.5 42812.5
 [811] 42813.5 42814.5 42815.5 42816.5 42817.5 42818.5 42819.5 42820.5 42821.5
 [820] 42822.5 42823.5 42824.5 42825.5 42826.5 42827.5 42828.5 42829.5 42830.5
 [829] 42831.5 42832.5 42833.5 42834.5 42835.5 42836.5 42837.5 42838.5 42839.5
 [838] 42840.5 42841.5 42842.5 42843.5 42844.5 42845.5 42846.5 42847.5 42848.5
 [847] 42849.5 42850.5 42851.5 42852.5 42853.5 42854.5 42855.5 42856.5 42857.5
 [856] 42858.5 42859.5 42860.5 42861.5 42862.5 42863.5 42864.5 42865.5 42866.5
 [865] 42867.5 42868.5 42869.5 42870.5 42871.5 42872.5 42873.5 42874.5 42875.5
 [874] 42876.5 42877.5 42878.5 42879.5 42880.5 42881.5 42882.5 42883.5 42884.5
 [883] 42885.5 42886.5 42887.5 42888.5 42889.5 42890.5 42891.5 42892.5 42893.5
 [892] 42894.5 42895.5 42896.5 42897.5 42898.5 42899.5 42900.5 42901.5 42902.5
 [901] 42903.5 42904.5 42905.5 42906.5 42907.5 42908.5 42909.5 42910.5 42911.5
 [910] 42912.5 42913.5 42914.5 42915.5 42916.5 42917.5 42918.5 42919.5 42920.5
 [919] 42921.5 42922.5 42923.5 42924.5 42925.5 42926.5 42927.5 42928.5 42929.5
 [928] 42930.5 42931.5 42932.5 42933.5 42934.5 42935.5 42936.5 42937.5 42938.5
 [937] 42939.5 42940.5 42941.5 42942.5 42943.5 42944.5 42945.5 42946.5 42947.5
 [946] 42948.5 42949.5 42950.5 42951.5 42952.5 42953.5 42954.5 42955.5 42956.5
 [955] 42957.5 42958.5 42959.5 42960.5 42961.5 42962.5 42963.5 42964.5 42965.5
 [964] 42966.5 42967.5 42968.5 42969.5 42970.5 42971.5 42972.5 42973.5 42974.5
 [973] 42975.5 42976.5 42977.5 42978.5 42979.5 42980.5 42981.5 42982.5 42983.5
 [982] 42984.5 42985.5 42986.5 42987.5 42988.5 42989.5 42990.5 42991.5 42992.5
 [991] 42993.5 42994.5 42995.5 42996.5 42997.5 42998.5 42999.5 43000.5 43001.5
[1000] 43002.5 43003.5 43004.5 43005.5 43006.5 43007.5 43008.5 43009.5 43010.5
[1009] 43011.5 43012.5 43013.5 43014.5 43015.5 43016.5 43017.5 43018.5 43019.5
[1018] 43020.5 43021.5 43022.5 43023.5 43024.5 43025.5 43026.5 43027.5 43028.5
[1027] 43029.5 43030.5 43031.5 43032.5 43033.5 43034.5 43035.5 43036.5 43037.5
[1036] 43038.5 43039.5 43040.5 43041.5 43042.5 43043.5 43044.5 43045.5 43046.5
[1045] 43047.5 43048.5 43049.5 43050.5 43051.5 43052.5 43053.5 43054.5 43055.5
[1054] 43056.5 43057.5 43058.5 43059.5 43060.5 43061.5 43062.5 43063.5 43064.5
[1063] 43065.5 43066.5 43067.5 43068.5 43069.5 43070.5 43071.5 43072.5 43073.5
[1072] 43074.5 43075.5 43076.5 43077.5 43078.5 43079.5 43080.5 43081.5 43082.5
[1081] 43083.5 43084.5 43085.5 43086.5 43087.5 43088.5 43089.5 43090.5 43091.5
[1090] 43092.5 43093.5 43094.5 43095.5 43096.5 43097.5 43098.5 43099.5 43100.5
[1099] 43101.5 43102.5 43103.5 43104.5 43105.5 43106.5 43107.5 43108.5 43109.5
[1108] 43110.5 43111.5 43112.5 43113.5 43114.5 43115.5 43116.5 43117.5 43118.5
[1117] 43119.5 43120.5 43121.5 43122.5 43123.5 43124.5 43125.5 43126.5 43127.5
[1126] 43128.5 43129.5 43130.5 43131.5 43132.5 43133.5 43134.5 43135.5 43136.5
[1135] 43137.5 43138.5 43139.5 43140.5 43141.5 43142.5 43143.5 43144.5 43145.5
[1144] 43146.5 43147.5 43148.5 43149.5 43150.5 43151.5 43152.5 43153.5 43154.5
[1153] 43155.5 43156.5 43157.5 43158.5 43159.5 43160.5 43161.5 43162.5 43163.5
[1162] 43164.5 43165.5 43166.5 43167.5 43168.5 43169.5 43170.5 43171.5 43172.5
[1171] 43173.5 43174.5 43175.5 43176.5 43177.5 43178.5 43179.5 43180.5 43181.5
[1180] 43182.5 43183.5 43184.5 43185.5 43186.5 43187.5 43188.5 43189.5 43190.5
[1189] 43191.5 43192.5 43193.5 43194.5 43195.5 43196.5 43197.5 43198.5 43199.5
[1198] 43200.5 43201.5 43202.5 43203.5 43204.5 43205.5 43206.5 43207.5 43208.5
[1207] 43209.5 43210.5 43211.5 43212.5 43213.5 43214.5 43215.5 43216.5 43217.5
[1216] 43218.5 43219.5 43220.5 43221.5 43222.5 43223.5 43224.5 43225.5 43226.5
[1225] 43227.5 43228.5 43229.5 43230.5 43231.5 43232.5 43233.5 43234.5 43235.5
[1234] 43236.5 43237.5 43238.5 43239.5 43240.5 43241.5 43242.5 43243.5 43244.5
[1243] 43245.5 43246.5 43247.5 43248.5 43249.5 43250.5 43251.5 43252.5 43253.5
[1252] 43254.5 43255.5 43256.5 43257.5 43258.5 43259.5 43260.5 43261.5 43262.5
[1261] 43263.5 43264.5 43265.5 43266.5 43267.5 43268.5 43269.5 43270.5 43271.5
[1270] 43272.5 43273.5 43274.5 43275.5 43276.5 43277.5 43278.5 43279.5 43280.5
[1279] 43281.5 43282.5 43283.5 43284.5 43285.5 43286.5 43287.5 43288.5 43289.5
[1288] 43290.5 43291.5 43292.5 43293.5 43294.5 43295.5 43296.5 43297.5 43298.5
[1297] 43299.5 43300.5 43301.5 43302.5 43303.5 43304.5 43305.5 43306.5 43307.5
[1306] 43308.5 43309.5 43310.5 43311.5 43312.5 43313.5 43314.5 43315.5 43316.5
[1315] 43317.5 43318.5 43319.5 43320.5 43321.5 43322.5 43323.5 43324.5 43325.5
[1324] 43326.5 43327.5 43328.5 43329.5 43330.5 43331.5 43332.5 43333.5 43334.5
[1333] 43335.5 43336.5 43337.5 43338.5 43339.5 43340.5 43341.5 43342.5 43343.5
[1342] 43344.5 43345.5 43346.5 43347.5 43348.5 43349.5 43350.5 43351.5 43352.5
[1351] 43353.5 43354.5 43355.5 43356.5 43357.5 43358.5 43359.5 43360.5 43361.5
[1360] 43362.5 43363.5 43364.5 43365.5 43366.5 43367.5 43368.5 43369.5 43370.5
[1369] 43371.5 43372.5 43373.5 43374.5 43375.5 43376.5 43377.5 43378.5 43379.5
[1378] 43380.5 43381.5 43382.5 43383.5 43384.5 43385.5 43386.5 43387.5 43388.5
[1387] 43389.5 43390.5 43391.5 43392.5 43393.5 43394.5 43395.5 43396.5 43397.5
[1396] 43398.5 43399.5 43400.5 43401.5 43402.5 43403.5 43404.5 43405.5 43406.5
[1405] 43407.5 43408.5 43409.5 43410.5 43411.5 43412.5 43413.5 43414.5 43415.5
[1414] 43416.5 43417.5 43418.5 43419.5 43420.5 43421.5 43422.5 43423.5 43424.5
[1423] 43425.5 43426.5 43427.5 43428.5 43429.5 43430.5 43431.5 43432.5 43433.5
[1432] 43434.5 43435.5 43436.5 43437.5 43438.5 43439.5 43440.5 43441.5 43442.5
[1441] 43443.5 43444.5 43445.5 43446.5 43447.5 43448.5 43449.5 43450.5 43451.5
[1450] 43452.5 43453.5 43454.5 43455.5 43456.5 43457.5 43458.5 43459.5 43460.5
[1459] 43461.5 43462.5 43463.5 43464.5 43465.5 43466.5 43467.5 43468.5 43469.5
[1468] 43470.5 43471.5 43472.5 43473.5 43474.5 43475.5 43476.5 43477.5 43478.5
[1477] 43479.5 43480.5 43481.5 43482.5 43483.5 43484.5 43485.5 43486.5 43487.5
[1486] 43488.5 43489.5 43490.5 43491.5 43492.5 43493.5 43494.5 43495.5 43496.5
[1495] 43497.5 43498.5 43499.5 43500.5 43501.5 43502.5 43503.5 43504.5 43505.5
[1504] 43506.5 43507.5 43508.5 43509.5 43510.5 43511.5 43512.5 43513.5 43514.5
[1513] 43515.5 43516.5 43517.5 43518.5 43519.5 43520.5 43521.5 43522.5 43523.5
[1522] 43524.5 43525.5 43526.5 43527.5 43528.5 43529.5 43530.5 43531.5 43532.5
[1531] 43533.5 43534.5 43535.5 43536.5 43537.5 43538.5 43539.5 43540.5 43541.5
[1540] 43542.5 43543.5 43544.5 43545.5 43546.5 43547.5 43548.5 43549.5 43550.5
[1549] 43551.5 43552.5 43553.5 43554.5 43555.5 43556.5 43557.5 43558.5 43559.5
[1558] 43560.5 43561.5 43562.5 43563.5 43564.5 43565.5 43566.5 43567.5 43568.5
[1567] 43569.5 43570.5 43571.5 43572.5 43573.5 43574.5 43575.5 43576.5 43577.5
[1576] 43578.5 43579.5 43580.5 43581.5 43582.5 43583.5 43584.5 43585.5 43586.5
[1585] 43587.5 43588.5 43589.5 43590.5 43591.5 43592.5 43593.5 43594.5 43595.5
[1594] 43596.5 43597.5 43598.5 43599.5 43600.5 43601.5 43602.5 43603.5 43604.5
[1603] 43605.5 43606.5 43607.5 43608.5 43609.5 43610.5 43611.5 43612.5 43613.5
[1612] 43614.5 43615.5 43616.5 43617.5 43618.5 43619.5 43620.5 43621.5 43622.5
[1621] 43623.5 43624.5 43625.5 43626.5 43627.5 43628.5 43629.5 43630.5 43631.5
[1630] 43632.5 43633.5 43634.5 43635.5 43636.5 43637.5 43638.5 43639.5 43640.5
[1639] 43641.5 43642.5 43643.5 43644.5 43645.5 43646.5 43647.5 43648.5 43649.5
[1648] 43650.5 43651.5 43652.5 43653.5 43654.5 43655.5 43656.5 43657.5 43658.5
[1657] 43659.5 43660.5 43661.5 43662.5 43663.5 43664.5 43665.5 43666.5 43667.5
[1666] 43668.5 43669.5 43670.5 43671.5 43672.5 43673.5 43674.5 43675.5 43676.5
[1675] 43677.5 43678.5 43679.5 43680.5 43681.5 43682.5 43683.5 43684.5 43685.5
[1684] 43686.5 43687.5 43688.5 43689.5 43690.5 43691.5 43692.5 43693.5 43694.5
[1693] 43695.5 43696.5 43697.5 43698.5 43699.5 43700.5 43701.5 43702.5 43703.5
[1702] 43704.5 43705.5 43706.5 43707.5 43708.5 43709.5 43710.5 43711.5 43712.5
[1711] 43713.5 43714.5 43715.5 43716.5 43717.5 43718.5 43719.5 43720.5 43721.5
[1720] 43722.5 43723.5 43724.5 43725.5 43726.5 43727.5 43728.5 43729.5 43730.5
[1729] 43731.5 43732.5 43733.5 43734.5 43735.5 43736.5 43737.5 43738.5 43739.5
[1738] 43740.5 43741.5 43742.5 43743.5 43744.5 43745.5 43746.5 43747.5 43748.5
[1747] 43749.5 43750.5 43751.5 43752.5 43753.5 43754.5 43755.5 43756.5 43757.5
[1756] 43758.5 43759.5 43760.5 43761.5 43762.5 43763.5 43764.5 43765.5 43766.5
[1765] 43767.5 43768.5 43769.5 43770.5 43771.5 43772.5 43773.5 43774.5 43775.5
[1774] 43776.5 43777.5 43778.5 43779.5 43780.5 43781.5 43782.5 43783.5 43784.5
[1783] 43785.5 43786.5 43787.5 43788.5 43789.5 43790.5 43791.5 43792.5 43793.5
[1792] 43794.5 43795.5 43796.5 43797.5 43798.5 43799.5 43800.5 43801.5 43802.5
[1801] 43803.5 43804.5 43805.5 43806.5 43807.5 43808.5 43809.5 43810.5 43811.5
[1810] 43812.5 43813.5 43814.5 43815.5 43816.5 43817.5 43818.5 43819.5 43820.5
[1819] 43821.5 43822.5 43823.5 43824.5 43825.5 43826.5 43827.5 43828.5 43829.5
[1828] 43830.5 43831.5 43832.5 43833.5 43834.5 43835.5 43836.5 43837.5 43838.5
[1837] 43839.5 43840.5 43841.5 43842.5 43843.5 43844.5 43845.5 43846.5 43847.5
[1846] 43848.5 43849.5 43850.5 43851.5 43852.5 43853.5 43854.5 43855.5 43856.5
[1855] 43857.5 43858.5 43859.5 43860.5 43861.5 43862.5 43863.5 43864.5 43865.5
[1864] 43866.5 43867.5 43868.5 43869.5 43870.5 43871.5 43872.5 43873.5 43874.5
[1873] 43875.5 43876.5 43877.5 43878.5 43879.5 43880.5 43881.5 43882.5 43883.5
[1882] 43884.5 43885.5 43886.5 43887.5 43888.5 43889.5 43890.5 43891.5 43892.5
[1891] 43893.5 43894.5 43895.5 43896.5 43897.5 43898.5 43899.5 43900.5 43901.5
[1900] 43902.5 43903.5 43904.5 43905.5 43906.5 43907.5 43908.5 43909.5 43910.5
[1909] 43911.5 43912.5 43913.5 43914.5 43915.5 43916.5 43917.5 43918.5 43919.5
[1918] 43920.5 43921.5 43922.5 43923.5 43924.5 43925.5 43926.5 43927.5 43928.5
[1927] 43929.5 43930.5 43931.5 43932.5 43933.5 43934.5 43935.5 43936.5 43937.5
[1936] 43938.5 43939.5 43940.5 43941.5 43942.5 43943.5 43944.5 43945.5 43946.5
[1945] 43947.5 43948.5 43949.5 43950.5 43951.5 43952.5 43953.5 43954.5 43955.5
[1954] 43956.5 43957.5 43958.5 43959.5 43960.5 43961.5 43962.5 43963.5 43964.5
[1963] 43965.5 43966.5 43967.5 43968.5 43969.5 43970.5 43971.5 43972.5 43973.5
[1972] 43974.5 43975.5 43976.5 43977.5 43978.5 43979.5 43980.5 43981.5 43982.5
[1981] 43983.5 43984.5 43985.5 43986.5 43987.5 43988.5 43989.5 43990.5 43991.5
[1990] 43992.5 43993.5 43994.5 43995.5 43996.5 43997.5 43998.5 43999.5 44000.5
[1999] 44001.5 44002.5 44003.5 44004.5 44005.5 44006.5 44007.5 44008.5 44009.5
[2008] 44010.5 44011.5 44012.5 44013.5 44014.5 44015.5 44016.5 44017.5 44018.5
[2017] 44019.5 44020.5 44021.5 44022.5 44023.5 44024.5 44025.5 44026.5 44027.5
[2026] 44028.5 44029.5 44030.5 44031.5 44032.5 44033.5 44034.5 44035.5 44036.5
[2035] 44037.5 44038.5 44039.5 44040.5 44041.5 44042.5 44043.5 44044.5 44045.5
[2044] 44046.5 44047.5 44048.5 44049.5 44050.5 44051.5 44052.5 44053.5 44054.5
[2053] 44055.5 44056.5 44057.5 44058.5 44059.5 44060.5 44061.5 44062.5 44063.5
[2062] 44064.5 44065.5 44066.5 44067.5 44068.5 44069.5 44070.5 44071.5 44072.5
[2071] 44073.5 44074.5 44075.5 44076.5 44077.5 44078.5 44079.5 44080.5 44081.5
[2080] 44082.5 44083.5 44084.5 44085.5 44086.5 44087.5 44088.5 44089.5 44090.5
[2089] 44091.5 44092.5 44093.5 44094.5 44095.5 44096.5 44097.5 44098.5 44099.5
[2098] 44100.5 44101.5 44102.5 44103.5 44104.5 44105.5 44106.5 44107.5 44108.5
[2107] 44109.5 44110.5 44111.5 44112.5 44113.5 44114.5 44115.5 44116.5 44117.5
[2116] 44118.5 44119.5 44120.5 44121.5 44122.5 44123.5 44124.5 44125.5 44126.5
[2125] 44127.5 44128.5 44129.5 44130.5 44131.5 44132.5 44133.5 44134.5 44135.5
[2134] 44136.5 44137.5 44138.5 44139.5 44140.5 44141.5 44142.5 44143.5 44144.5
[2143] 44145.5 44146.5 44147.5 44148.5 44149.5 44150.5 44151.5 44152.5 44153.5
[2152] 44154.5 44155.5 44156.5 44157.5 44158.5 44159.5 44160.5 44161.5 44162.5
[2161] 44163.5 44164.5 44165.5 44166.5 44167.5 44168.5 44169.5 44170.5 44171.5
[2170] 44172.5 44173.5 44174.5 44175.5 44176.5 44177.5 44178.5 44179.5 44180.5
[2179] 44181.5 44182.5 44183.5 44184.5 44185.5 44186.5 44187.5 44188.5 44189.5
[2188] 44190.5 44191.5 44192.5 44193.5 44194.5 44195.5 44196.5 44197.5 44198.5
[2197] 44199.5 44200.5 44201.5 44202.5 44203.5 44204.5 44205.5 44206.5 44207.5
[2206] 44208.5 44209.5 44210.5 44211.5 44212.5 44213.5 44214.5 44215.5 44216.5
[2215] 44217.5 44218.5 44219.5 44220.5 44221.5 44222.5 44223.5 44224.5 44225.5
[2224] 44226.5 44227.5 44228.5 44229.5 44230.5 44231.5 44232.5 44233.5 44234.5
[2233] 44235.5 44236.5 44237.5 44238.5 44239.5 44240.5 44241.5 44242.5 44243.5
[2242] 44244.5 44245.5 44246.5 44247.5 44248.5 44249.5 44250.5 44251.5 44252.5
[2251] 44253.5 44254.5 44255.5 44256.5 44257.5 44258.5 44259.5 44260.5 44261.5
[2260] 44262.5 44263.5 44264.5 44265.5 44266.5 44267.5 44268.5 44269.5 44270.5
[2269] 44271.5 44272.5 44273.5 44274.5 44275.5 44276.5 44277.5 44278.5 44279.5
[2278] 44280.5 44281.5 44282.5 44283.5 44284.5 44285.5 44286.5 44287.5 44288.5
[2287] 44289.5 44290.5 44291.5 44292.5 44293.5 44294.5 44295.5 44296.5 44297.5
[2296] 44298.5 44299.5 44300.5 44301.5 44302.5 44303.5 44304.5 44305.5 44306.5
[2305] 44307.5 44308.5 44309.5 44310.5 44311.5 44312.5 44313.5 44314.5 44315.5
[2314] 44316.5 44317.5 44318.5 44319.5 44320.5 44321.5 44322.5 44323.5 44324.5
[2323] 44325.5 44326.5 44327.5 44328.5 44329.5 44330.5 44331.5 44332.5 44333.5
[2332] 44334.5 44335.5 44336.5 44337.5 44338.5 44339.5 44340.5 44341.5 44342.5
[2341] 44343.5 44344.5 44345.5 44346.5 44347.5 44348.5 44349.5 44350.5 44351.5
[2350] 44352.5 44353.5 44354.5 44355.5 44356.5 44357.5 44358.5 44359.5 44360.5
[2359] 44361.5 44362.5 44363.5 44364.5 44365.5 44366.5 44367.5 44368.5 44369.5
[2368] 44370.5 44371.5 44372.5 44373.5 44374.5 44375.5 44376.5 44377.5 44378.5
[2377] 44379.5 44380.5 44381.5 44382.5 44383.5 44384.5 44385.5 44386.5 44387.5
[2386] 44388.5 44389.5 44390.5 44391.5 44392.5 44393.5 44394.5 44395.5 44396.5
[2395] 44397.5 44398.5 44399.5 44400.5 44401.5 44402.5 44403.5 44404.5 44405.5
[2404] 44406.5 44407.5 44408.5 44409.5 44410.5 44411.5 44412.5 44413.5 44414.5
[2413] 44415.5 44416.5 44417.5 44418.5 44419.5 44420.5 44421.5 44422.5 44423.5
[2422] 44424.5 44425.5 44426.5 44427.5 44428.5 44429.5 44430.5 44431.5 44432.5
[2431] 44433.5 44434.5 44435.5 44436.5 44437.5 44438.5 44439.5 44440.5 44441.5
[2440] 44442.5 44443.5 44444.5 44445.5 44446.5 44447.5 44448.5 44449.5 44450.5
[2449] 44451.5 44452.5 44453.5 44454.5 44455.5 44456.5 44457.5 44458.5 44459.5
[2458] 44460.5 44461.5 44462.5 44463.5 44464.5 44465.5 44466.5 44467.5 44468.5
[2467] 44469.5 44470.5 44471.5 44472.5 44473.5 44474.5 44475.5 44476.5 44477.5
[2476] 44478.5 44479.5 44480.5 44481.5 44482.5 44483.5 44484.5 44485.5 44486.5
[2485] 44487.5 44488.5 44489.5 44490.5 44491.5 44492.5 44493.5 44494.5 44495.5
[2494] 44496.5 44497.5 44498.5 44499.5 44500.5 44501.5 44502.5 44503.5 44504.5
[2503] 44505.5 44506.5 44507.5 44508.5 44509.5 44510.5 44511.5 44512.5 44513.5
[2512] 44514.5 44515.5 44516.5 44517.5 44518.5 44519.5 44520.5 44521.5 44522.5
[2521] 44523.5 44524.5 44525.5 44526.5 44527.5 44528.5 44529.5 44530.5 44531.5
[2530] 44532.5 44533.5 44534.5 44535.5 44536.5 44537.5 44538.5 44539.5 44540.5
[2539] 44541.5 44542.5 44543.5 44544.5 44545.5 44546.5 44547.5 44548.5 44549.5
[2548] 44550.5 44551.5 44552.5 44553.5 44554.5 44555.5 44556.5 44557.5 44558.5
[2557] 44559.5
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 2557
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get oxygen
dname <- "o2b"

oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
[1]  336  289 2557
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA

dim(oxy_array)
[1]  336  289 2557
str(dim(oxy_array))
 int [1:3] 336 289 2557
# The third slot is the date index

# Loop through all "dates", put into a list 
dlist <- list()

for(i in 1:length(months)) {
  
  oxy_sub <- oxy_array[, , i]
    
  dlist[[i]] <- oxy_sub
  
}

# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)
List of 2557
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 1_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 2_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 3_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
 $ 4_2015 : num [1:336, 1:289] NA NA NA NA NA NA NA NA NA NA ...
  [list output truncated]
# Create data holding object
oxy_data_list <- list()

# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat_1$month_year)) {
  
  # Set plot limits
  ymin = 54; ymax = 58; xmin = 12; xmax = 22

  # Subset a month-year combination
  oxy_slice <- dlist[[i]]
  
  # Create raster for that year (i)
  r <- raster(t(oxy_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
              crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r <- flip(r, direction = 'y')
  
  plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))

  # Filter the same year (i) in the cpue data and select only coordinates
  d_slice <- dat_1 %>% filter(month_year == i) %>% dplyr::select(lon, lat)
  
  # Make into a SpatialPoints object
  data_sp <- SpatialPoints(d_slice)
  
  # Extract raster value (oxygen)
  rasValue <- raster::extract(r, data_sp)
  
  # Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for plot)
  df <- as.data.frame(data_sp)
  
  # Add in the raster value in the df holding the coordinates for the cpue data
  d_slice$oxy <- rasValue
  
  # Add in which year
  d_slice$month_year <- i

  # Now the unit of oxygen is mmol/m3. I want it to be ml/L. The original model is in unit ml/L
  # and it's been converted by the data host. Since it was converted without accounting for
  # pressure or temperature, I can simply use the following conversion factor:
  # 1 ml/l = 103/22.391 = 44.661 μmol/l -> 1 ml/l = 0.044661 mmol/l = 44.661 mmol/m^3 -> 0.0223909 ml/l = 1mmol/m^3
  # https://ocean.ices.dk/tools/unitconversion.aspx

  d_slice$oxy <- d_slice$oxy * 0.0223909
    
  # Add each years' data in the list
  oxy_data_list[[i]] <- d_slice

}
filter: removed 289 rows (95%), 16 rows remaining

filter: removed 275 rows (90%), 30 rows remaining

filter: removed 286 rows (94%), 19 rows remaining

filter: removed 269 rows (88%), 36 rows remaining

filter: removed 288 rows (94%), 17 rows remaining

filter: removed 274 rows (90%), 31 rows remaining

filter: removed 299 rows (98%), 6 rows remaining

filter: removed 280 rows (92%), 25 rows remaining

filter: removed 265 rows (87%), 40 rows remaining

filter: removed 288 rows (94%), 17 rows remaining

filter: removed 264 rows (87%), 41 rows remaining

filter: removed 278 rows (91%), 27 rows remaining

# Now create a data frame from the list of all annual values
big_dat_oxy <- dplyr::bind_rows(oxy_data_list)

sort(unique(big_dat_oxy$month_year))
 [1] "11_2015" "11_2016" "11_2017" "11_2018" "11_2019" "11_2020" "11_2021"
 [8] "2_2016"  "2_2017"  "2_2018"  "2_2020"  "2_2021" 

New oxygen

# 2022 Q1
# Note this is from the projection, not the re-analysis because it hasn't been updated
ncin <- nc_open(paste0(home, "/data/NEMO/BAL-ERGOM_BGC-MonthlyMeans-202202.nc"))

print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/NEMO/BAL-ERGOM_BGC-MonthlyMeans-202202.nc (NC_FORMAT_NETCDF4_CLASSIC):

     13 variables (excluding dimension variables):
        float chl[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mg m-3
            long_name: Chlorophyll_a_Concentration
            standard_name: mass_concentration_of_chlorophyll_a_in_sea_water
            unit_long: mass per unit volume
            missing_value: -999
        float phyc[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Phytoplankton Concentration
            standard_name: mole_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float o2b[lon,lat,time]   (Chunking: [763,774,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Dissolved_Oxygen_Concentration
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            unit_long: millimole O2 per unit volume
            missing_value: -999
        float no3[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Nitrate_Concentration
            standard_name: mole_concentration_of_nitrate_in_sea_water
            unit_long: millimole N per unit volume
            missing_value: -999
        float kd[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: 1/m
            long_name: Light Attenuation Coefficient
            standard_name: volume_attenuation_coefficient_of_downwelling_radiative_flux_in_sea_water
            unit_long: 1/m
            missing_value: -999
        float o2[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Dissolved_Oxygen_Concentration
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            unit_long: millimole O2 per unit volume
            missing_value: -999
        float pH[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: 1
            long_name: sea_water_ph
            standard_name: sea_water_ph_reported_on_total_scale
            unit_long: 1
            missing_value: -999
        float po4[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Phosphate_Concentration
            standard_name: mole_concentration_of_phosphate_in_sea_water
            unit_long: millimole P per unit volume
            missing_value: -999
        float zooc[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Zooplankton Concentration
            standard_name: mole_concentration_of_zooplankton_expressed_as_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float spCO2[lon,lat,time]   (Chunking: [763,774,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: Pa
            long_name: surface_partial_pressure_of_carbon_dioxide
            standard_name: surface_partial_pressure_of_carbon_dioxide_in_sea_water
            unit_long: pascal
            missing_value: -999
        float nh4[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Ammonium_Concentration
            standard_name: mole_concentration_of_ammonium_in_sea_water
            unit_long: millimole N per unit volume
            missing_value: -999
        float dissic[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mol m-3
            long_name: Dissolved_Inorganic_Carbon_Concentration
            standard_name: mole_concentration_of_dissolved_inorganic_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float si[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Silica_Concentration
            standard_name: mole_concentration_of_silicate_in_sea_water
            unit_long: millimole Silicate per unit volume
            missing_value: -999

     4 dimensions:
        time  Size:1   *** is unlimited *** 
            units: days since 1900-01-01 00:00:00
            long_name: time
            standard_name: time
            axis: T
        depth  Size:56 
            units: m
            comment: 
            unit_long: meters
            reference: sea_level SeaDataNet L111
            positive: down
            uncertainty: 
            valid_min: 0
            QC_procedure: 1
            long_name: Depth of each measurement
            standard_name: depth
            QC_indicator: 1
            valid_max: 12000
            axis: Z
        lon  Size:763 
            units: degrees_east
            long_name: Longitude of each location
            standard_name: longitude
            unit_long: degrees East
            axis: X
        lat  Size:774 
            units: degrees_north
            long_name: Latitude of each location
            standard_name: latitude
            unit_long: degrees North
            axis: Y

    18 global attributes:
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        title: CMEMS ERGOM monthly integrated model fields
        Conventions: CF-1.0
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        institution: Baltic MFC, PU Swedish Meteorological and Hydrological Institute
        compression: yes
        stop_date: 2022-02-15 03:29:27
        creation_date: 2023-09-30 01:48:04
        netcdf_version_id: 4.4.1.1 of Oct 24 2018 14:14:02 $
        easternmost_longitude: 30.2086563110352
        northernmost_latitude: 65.8909912109375
        westernmost_longitude: 9.04158210754395
        southernmost_latitude: 53.0082969665527
        start_date: 2022-02-15 03:29:27
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 9.041582 9.069360 9.097137 9.124915 9.152693 9.180470
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163
# Get time
time <- ncvar_get(ncin,"time")
time
[1] 44605.15
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 1
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get oxygen
dname <- "o2b"

oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
[1] 763 774
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA

# Create raster for that year (i)
r <- raster(t(oxy_array), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
            crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
# Flip...
r <- flip(r, direction = 'y')
  
plot(r, main = paste("2022_q1"), ylim = c(ymin, ymax), xlim = c(xmin, xmax))

# Filter the same year (i) in the cpue data and select only coordinates
dat_2_q1 <- dat_2 %>% dplyr::select(lon, lat)
  
# Make into a SpatialPoints object
dat_2_q1_sp <- SpatialPoints(dat_2_q1)
  
# Extract raster value (oxygen)
rasValue <- raster::extract(r, dat_2_q1_sp)
  
dat_2_q1$oxy <- rasValue
  
# Add in which year
dat_2_q1$month_year <- "2_2022"

head(dat_2_q1)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.4  55.7  307. 2_2022    
2  14.4  55.7  287. 2_2022    
3  14.7  55.7  282. 2_2022    
4  15.9  55.8  242. 2_2022    
5  13.7  55.2  361. 2_2022    
6  14.7  55.6  240. 2_2022    
dat_2_q1$oxy <- dat_2_q1$oxy * 0.0223909
    
hist(dat_2_q1$oxy)

# Now do Q4
ncin <- nc_open(paste0(home, "/data/NEMO/BAL-ERGOM_BGC-MonthlyMeans-202211.nc"))

print(ncin)
File /Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/NEMO/BAL-ERGOM_BGC-MonthlyMeans-202211.nc (NC_FORMAT_NETCDF4_CLASSIC):

     13 variables (excluding dimension variables):
        float chl[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mg m-3
            long_name: Chlorophyll_a_Concentration
            standard_name: mass_concentration_of_chlorophyll_a_in_sea_water
            unit_long: mass per unit volume
            missing_value: -999
        float phyc[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Phytoplankton Concentration
            standard_name: mole_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float o2b[lon,lat,time]   (Chunking: [763,774,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Dissolved_Oxygen_Concentration
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            unit_long: millimole O2 per unit volume
            missing_value: -999
        float no3[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Nitrate_Concentration
            standard_name: mole_concentration_of_nitrate_in_sea_water
            unit_long: millimole N per unit volume
            missing_value: -999
        float kd[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: 1/m
            long_name: Light Attenuation Coefficient
            standard_name: volume_attenuation_coefficient_of_downwelling_radiative_flux_in_sea_water
            unit_long: 1/m
            missing_value: -999
        float o2[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Dissolved_Oxygen_Concentration
            standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
            unit_long: millimole O2 per unit volume
            missing_value: -999
        float pH[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: 1
            long_name: sea_water_ph
            standard_name: sea_water_ph_reported_on_total_scale
            unit_long: 1
            missing_value: -999
        float po4[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Phosphate_Concentration
            standard_name: mole_concentration_of_phosphate_in_sea_water
            unit_long: millimole P per unit volume
            missing_value: -999
        float zooc[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Zooplankton Concentration
            standard_name: mole_concentration_of_zooplankton_expressed_as_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float spCO2[lon,lat,time]   (Chunking: [763,774,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: Pa
            long_name: surface_partial_pressure_of_carbon_dioxide
            standard_name: surface_partial_pressure_of_carbon_dioxide_in_sea_water
            unit_long: pascal
            missing_value: -999
        float nh4[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Ammonium_Concentration
            standard_name: mole_concentration_of_ammonium_in_sea_water
            unit_long: millimole N per unit volume
            missing_value: -999
        float dissic[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mol m-3
            long_name: Dissolved_Inorganic_Carbon_Concentration
            standard_name: mole_concentration_of_dissolved_inorganic_carbon_in_sea_water
            unit_long: millimole C per unit volume
            missing_value: -999
        float si[lon,lat,depth,time]   (Chunking: [191,194,14,1])  (Compression: shuffle,level 1)
            _FillValue: -999
            units: mmol m-3
            long_name: Silica_Concentration
            standard_name: mole_concentration_of_silicate_in_sea_water
            unit_long: millimole Silicate per unit volume
            missing_value: -999

     4 dimensions:
        time  Size:1   *** is unlimited *** 
            units: days since 1900-01-01 00:00:00
            long_name: time
            standard_name: time
            axis: T
        depth  Size:56 
            units: m
            comment: 
            unit_long: meters
            reference: sea_level SeaDataNet L111
            positive: down
            uncertainty: 
            valid_min: 0
            QC_procedure: 1
            long_name: Depth of each measurement
            standard_name: depth
            QC_indicator: 1
            valid_max: 12000
            axis: Z
        lon  Size:763 
            units: degrees_east
            long_name: Longitude of each location
            standard_name: longitude
            unit_long: degrees East
            axis: X
        lat  Size:774 
            units: degrees_north
            long_name: Latitude of each location
            standard_name: latitude
            unit_long: degrees North
            axis: Y

    18 global attributes:
        comment: Data on cropped native product grid. Horizontal velocities destagged
        grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
        title: CMEMS ERGOM monthly integrated model fields
        Conventions: CF-1.0
        source: CMEMS BAL MFC NEMO model output converted to NetCDF
        contact: servicedesk.cmems@mercator-ocean.eu
        references: https://marine.copernicus.eu/
        file_quality_index: 1
        institution: Baltic MFC, PU Swedish Meteorological and Hydrological Institute
        compression: yes
        stop_date: 2022-11-27 06:00:00
        creation_date: 2023-11-08 23:27:07
        netcdf_version_id: 4.4.1.1 of Oct 24 2018 14:14:02 $
        easternmost_longitude: 30.2086563110352
        northernmost_latitude: 65.8909912109375
        westernmost_longitude: 9.04158210754395
        southernmost_latitude: 53.0082969665527
        start_date: 2022-11-27 06:00:00
# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)
[1] 9.041582 9.069360 9.097137 9.124915 9.152693 9.180470
lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)
[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163
# Get time
time <- ncvar_get(ncin,"time")
time
[1] 44890.25
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
[1] 1
tunits
$hasatt
[1] TRUE

$value
[1] "days since 1900-01-01 00:00:00"
# Get oxygen
dname <- "o2b"

oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
[1] 763 774
# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA

# Create raster for that year (i)
r <- raster(t(oxy_array), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
            crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
# Flip...
r <- flip(r, direction = 'y')
  
plot(r, main = paste("2022_q4"), ylim = c(ymin, ymax), xlim = c(xmin, xmax))

# Filter the same year (i) in the cpue data and select only coordinates
dat_2_q4 <- dat_2 %>% dplyr::select(lon, lat)
  
# Make into a SpatialPoints object
dat_2_q4_sp <- SpatialPoints(dat_2_q4)
  
# Extract raster value (oxygen)
rasValue <- raster::extract(r, dat_2_q4_sp)
  
dat_2_q4$oxy <- rasValue
  
# Add in which year
dat_2_q4$month_year <- "11_2022"

head(dat_2_q4)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.4  55.7  328. 11_2022   
2  14.4  55.7  324. 11_2022   
3  14.7  55.7  207. 11_2022   
4  15.9  55.8  167. 11_2022   
5  13.7  55.2  258. 11_2022   
6  14.7  55.6  161. 11_2022   
dat_2_q4$oxy <- dat_2_q4$oxy * 0.0223909
    
hist(dat_2_q4$oxy)

New oxygen

head(dat_2_q1)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.4  55.7  6.88 2_2022    
2  14.4  55.7  6.43 2_2022    
3  14.7  55.7  6.31 2_2022    
4  15.9  55.8  5.42 2_2022    
5  13.7  55.2  8.08 2_2022    
6  14.7  55.6  5.36 2_2022    
head(dat_2_q4)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.4  55.7  7.35 11_2022   
2  14.4  55.7  7.25 11_2022   
3  14.7  55.7  4.64 11_2022   
4  15.9  55.8  3.74 11_2022   
5  13.7  55.2  5.78 11_2022   
6  14.7  55.6  3.60 11_2022   
hist(dat_2_q1$oxy)

hist(dat_2_q4$oxy)

Join with stomach data

d_wide4 <- d_wide3

d_wide4 <- d_wide4 %>% 
  mutate(month_2 = ifelse(quarter == 1, 2, 11), # most common months
         month_year = paste(month_2, year, sep = "_")) %>% 
  dplyr::select(-month_2)
mutate: new variable 'month_2' (double) with 2 unique values and 0% NA
        new variable 'month_year' (character) with 14 unique values and 0% NA
unique(d_wide4$month_year)
 [1] "11_2015" "2_2016"  "11_2016" "2_2017"  "11_2017" "2_2018"  "11_2018"
 [8] "11_2019" "2_2020"  "11_2020" "2_2021"  "11_2021" "2_2022"  "11_2022"
unique(big_dat_oxy$month_year)
 [1] "11_2015" "2_2016"  "11_2016" "2_2017"  "11_2017" "2_2018"  "11_2018"
 [8] "11_2019" "2_2020"  "11_2020" "2_2021"  "11_2021"
# Add saduria
d_wide4 <- d_wide4 %>% left_join(dat %>% dplyr::select(X, Y, density_saduria), by = c("X", "Y"))
left_join: added one column (density_saduria)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     25,830    (includes duplicates)
           >                 ========
           > rows total       25,830
# Add oxygen (bind_rows first)
head(big_dat_oxy)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.6  55.6  3.97 11_2015   
2  14.5  55.7  5.44 11_2015   
3  14.4  55.7  6.25 11_2015   
4  17.4  56.1  5.59 11_2015   
5  18.4  56.3  7.01 11_2015   
6  18.6  56.4  4.83 11_2015   
head(dat_2_q1)
# A tibble: 6 × 4
    lon   lat   oxy month_year
  <dbl> <dbl> <dbl> <chr>     
1  14.4  55.7  6.88 2_2022    
2  14.4  55.7  6.43 2_2022    
3  14.7  55.7  6.31 2_2022    
4  15.9  55.8  5.42 2_2022    
5  13.7  55.2  8.08 2_2022    
6  14.7  55.6  5.36 2_2022    
oxy_all <- bind_rows(big_dat_oxy, dat_2_q1, dat_2_q4)

ggplot(oxy_all, aes(lon, lat, color = oxy)) + 
  geom_point() + 
  facet_wrap(~month_year)

d_wide4 <- d_wide4 %>% left_join(oxy_all, by = c("lon", "lat", "month_year"))
left_join: added one column (oxy)
           > rows only in x        0
           > rows only in y  (    43)
           > matched rows     30,362    (includes duplicates)
           >                 ========
           > rows total       30,362
ggplot(d_wide4, aes(lon, lat, color = oxy)) + 
  geom_point() + 
  facet_wrap(~month_year)

d_wide4 <- d_wide4 %>% dplyr::select(-month_year)

Save

write_csv(d_wide4, paste0(home, "/data/clean/aggregated_stomach_data.csv"))